NATRE microstructure dataset
Contents
NATRE microstructure dataset#
Study NATRE dataset, reproduce Ferrari & Polzin (2005)
%load_ext watermark
%matplotlib inline
import glob
import cf_xarray as cfxr
import dcpy
import distributed
import gsw_xarray
import hvplot.xarray
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import tqdm
import xfilter
from IPython.display import Image
import eddydiff as ed
import xarray as xr
xr.set_options(keep_attrs=True)
plt.rcParams["figure.dpi"] = 140
plt.rcParams["savefig.dpi"] = 200
plt.style.use("ggplot")
%watermark -iv
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
from pandas import Int64Index as NumericIndex
numpy : 1.22.3
tqdm : 4.64.0
distributed: 2022.4.0
hvplot : 0.7.3
cf_xarray : 0.7.1.dev16+g7305fd8
gsw_xarray : 0.3.0
xarray : 2022.6.0
matplotlib : 3.5.1
eddydiff : 0.1
dcpy : 0.2.1.dev11+gebbbf87.d20220428
xfilter : 0.2.1.dev6+g5fdc006
if "client" in locals():
client.cluster.close()
client.close()
client = distributed.Client(n_workers=2, threads_per_worker=2)
client
Client
Client-8948807e-35dd-11ed-95d3-3af9d394f1c6
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: http://127.0.0.1:8787/status |
Cluster Info
LocalCluster
f257c023
| Dashboard: http://127.0.0.1:8787/status | Workers: 2 |
| Total threads: 4 | Total memory: 16.00 GiB |
| Status: running | Using processes: True |
Scheduler Info
Scheduler
Scheduler-35410ff6-70d8-4bb7-85bd-0a821e8abd2a
| Comm: tcp://127.0.0.1:54020 | Workers: 2 |
| Dashboard: http://127.0.0.1:8787/status | Total threads: 4 |
| Started: Just now | Total memory: 16.00 GiB |
Workers
Worker: 0
| Comm: tcp://127.0.0.1:54032 | Total threads: 2 |
| Dashboard: http://127.0.0.1:54034/status | Memory: 8.00 GiB |
| Nanny: tcp://127.0.0.1:54024 | |
| Local directory: /Users/dcherian/work/eddydiff/notebooks/dask-worker-space/worker-c83c9tud | |
Worker: 1
| Comm: tcp://127.0.0.1:54031 | Total threads: 2 |
| Dashboard: http://127.0.0.1:54033/status | Memory: 8.00 GiB |
| Nanny: tcp://127.0.0.1:54023 | |
| Local directory: /Users/dcherian/work/eddydiff/notebooks/dask-worker-space/worker-5g_85fsd | |
Image("../images/natre-large-scale.png")
Read dataset#
created using eddydiff.natre.combine_natre_files()
natre = ed.natre.read_natre(load=True)
13 40
0 values == 0
28090
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in get_gradient_sign
return self.ufunc(*args, **kwargs)
ed.plot.histogram_turb_estimates(natre)
\(R_ρ, K_T, K_ρ\)#
(
natre.Rρ.sel(pres=slice(250, None))
.where(np.abs(natre.Rρ) < 100)
.mean(["latitude", "longitude"])
.reset_coords(drop=True)
.hvplot()
).opts(show_grid=True)
ds = xr.Dataset(
{"KrhoKt": (natre.Krho / natre.Kt), "Rρ": natre.Rρ.where(np.abs(natre.Rρ) < 5)}
).sel(pres=slice(250, None))
ds.plot.scatter(y="KrhoKt", x="Rρ")
plt.yscale("log")
import colorcet
import datashader as dshdr
from holoviews.operation.datashader import datashade
df = ds.to_dataframe()
# df["KrhoKt"] = np.log10(df.KrhoKt)
pts = hv.Points(df, ["KrhoKt", "Rρ"])
datashade(pts, cmap=colorcet.fire).opts(width=600, height=600)
TS plot#
hdl, axes = dcpy.oceans.TSplot(
natre.salt[::2],
natre.theta[::2],
Pref=1000,
# rho_levels=bins,
Sbins=np.arange(35, 35.5, 0.025),
Tbins=np.arange(4, 10, 0.3),
hexbin=False,
plot_kwargs={"alpha": 0.01},
)
(
natre["chi"].rolling(depth=100, center=True, min_periods=1).mean().compute()
).plot.line(
y="depth",
col="longitude",
row="latitude",
yincrease=False,
xlim=[1e-12, 1e-7],
ylim=(2000, 0),
xscale="log",
)
tornado.application - ERROR - Exception in callback functools.partial(<bound method IOLoop._discard_future_result of <zmq.eventloop.ioloop.ZMQIOLoop object at 0x7f70c8187340>>, <Task finished name='Task-481' coro=<Cluster._sync_cluster_info() done, defined at /home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py:104> exception=OSError('Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s')>)
Traceback (most recent call last):
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 284, in connect
comm = await asyncio.wait_for(
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/asyncio/tasks.py", line 498, in wait_for
raise exceptions.TimeoutError()
asyncio.exceptions.TimeoutError
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 741, in _run_callback
ret = callback()
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 765, in _discard_future_result
future.result()
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py", line 105, in _sync_cluster_info
await self.scheduler_comm.set_metadata(
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 785, in send_recv_from_rpc
comm = await self.live_comm()
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 742, in live_comm
comm = await connect(
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 308, in connect
raise OSError(
OSError: Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s
tornado.application - ERROR - Exception in callback functools.partial(<bound method IOLoop._discard_future_result of <zmq.eventloop.ioloop.ZMQIOLoop object at 0x7f70c8187340>>, <Task finished name='Task-488' coro=<Cluster._sync_cluster_info() done, defined at /home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py:104> exception=OSError('Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s')>)
Traceback (most recent call last):
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 284, in connect
comm = await asyncio.wait_for(
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/asyncio/tasks.py", line 498, in wait_for
raise exceptions.TimeoutError()
asyncio.exceptions.TimeoutError
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 741, in _run_callback
ret = callback()
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 765, in _discard_future_result
future.result()
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py", line 105, in _sync_cluster_info
await self.scheduler_comm.set_metadata(
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 785, in send_recv_from_rpc
comm = await self.live_comm()
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 742, in live_comm
comm = await connect(
File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 308, in connect
raise OSError(
OSError: Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s
<xarray.plot.facetgrid.FacetGrid at 0x7f70a01c4e50>
Procedure#
The large-scale average operator ⟨⟩ represents
a horizontal average over the survey lateral scale,
a vertical average over O(100) m, and
a time average over the 18-day survey.
The mean fields are derived by averaging all variables along neutral-density surfaces \(γ_n\)
Mean vertical gradients \(∂_z θ_m\), \(∂_z S_m\), and \(∂_z b_m\), are calculated from O(100)-m linear fits to \(θ_m(z_n)\), \(S_m(z_n)\), and \(b_m(z_n)\), where \(z_n\) is the mean depth of each surface \(γ_n\).
large-scale temperature gradient \(|∇_n θ_m|\) is estimated by fitting a plane to the 100 stations in the 400 km × 400 km large-scale survey grid and rms is computed as the standard deviation of the departures from the plane fit.
In theory, these bins are approx. 100m apart in neutral density
Choose bins#
# bins = (
# natre.gamma_n.mean(["latitude", "longitude"])
# .interp(depth=np.arange(150, 2001, 100))
# .dropna("depth")
# .data
# )
bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
natre.gamma_n.stack({"latlon": ["latitude", "longitude"]}).drop("latlon").cf.plot(
hue="latlon",
y="depth",
lw=0.5,
add_legend=False,
yincrease=False,
)
dcpy.plots.linex(bins)
Make estimate#
Vertical gradients: ed.sections.fit1D
Mean vertical gradients \(∂_z θ_m\), \(∂_z S_m\), and \(∂_z b_m\), are calculated from O(100)-m linear fits to \(θ_m(z_n)\), \(S_m(z_n)\), and \(b_m(z_n)\), where \(z_n\) is the mean depth of each surface \(γ_n\).
I don’t understand this. The bins are O(100m) apart, so how do you “fit” straight lines over O(100)-m to a profile that has points every O(100)m.
Horizontal gradients: ed.sections.fit2D
large-scale temperature gradient \(|∇_n θ_m|\) is estimated by fitting a plane to the 100 stations in the 400 km × 400 km large-scale survey grid and rms is computed as the standard deviation of the departures from the plane fit.
Ke is totally wrong!, I need to fix my plane fitting
Bin-average in density bins#
micro = ed.sections.bin_average_vertical(
natre.reset_coords("pres"), "neutral_density", bins
)
micro
<xarray.Dataset>
Dimensions: (gamma_n: 17, bounds: 2)
Coordinates:
* gamma_n (gamma_n) float64 26.79 26.96 27.1 ... 27.93 27.95 27.96
pres (gamma_n) float64 303.1 407.2 ... 1.826e+03 1.902e+03
reference_pressure int64 1000
num_obs (gamma_n) int64 19763 19424 17093 ... 18916 21881 10032
gamma_n_bounds (bounds, gamma_n) float64 26.69 26.88 ... 27.96 27.97
Dimensions without coordinates: bounds
Data variables: (12/18)
chi (gamma_n) float64 6.112e-09 3.64e-09 ... 2.707e-10
eps (gamma_n) float64 6.149e-10 5.53e-10 ... 1.16e-10
salt (gamma_n) float64 36.16 35.87 35.68 ... 35.13 35.1 35.09
temp (gamma_n) float64 15.64 13.83 12.45 ... 4.776 4.47 4.271
theta (gamma_n) float64 15.75 13.92 12.52 ... 4.393 4.187
pden (gamma_n) float64 1.031e+03 1.031e+03 ... 1.032e+03
... ...
Kt (gamma_n) float64 2.416e-05 2.335e-05 ... 1.511e-05
KtTz (gamma_n) float64 1.398e-07 1.077e-07 ... 3.357e-08
dTdz_m (gamma_n) float64 0.01861 0.01519 ... 0.002752 0.002534
N2_m (gamma_n) float64 1.95e-05 1.679e-05 ... 1.355e-06
Krho_m (gamma_n) float64 6.307e-06 6.589e-06 ... 1.713e-05
Kt_m (gamma_n) float64 8.828e-06 7.89e-06 ... 2.107e-05
Attributes: (12/13)
Conventions: CF-1.6
netcdf_version: 4
project: North Atlantic Tracer Release Experiment (NATRE)
expocode: 32OC250_4
cast_number: 3.0
title: Microstructure profiler data from the ship Oceanus...
... ...
latitude: 27.533166666666666
longitude: -30.723333333333333
chief_scientist: Raymond W. Schmitt
data_originator: Polzin
institution: WHOI
data_assembly_center: CCHDOReplicate some Ferrari & Polzin (2005) figures. \(K_ρ, ε\) look OK. \(K_e\) is wrong, something is wrong with my slope estimate / plane fitting
f, ax = plt.subplots(1, 3, sharey=True, constrained_layout=True)
micro.drop_vars("pres").Krho_m.cf.plot.step(xscale="log", xlim=(4e-6, 1e-4), ax=ax[0])
micro.drop_vars("pres").eps.cf.plot.step(ax=ax[1], xscale="log")
# chidens["Ke"].cf.plot(ax=ax[2])
[<matplotlib.lines.Line2D at 0x7f706d7c6e20>]
Constructing mean profiles#
I experimented with constructing a mean profile using a large number of γ_n bins, and then fitting straight lines in approx 100m bins.
I have instead decided to do a similar fitting procedure in ed.sections.fit1D for each chosen O(100m) wide density bin. Otherwise I would need to figure out how to go from \(θ_z\) in this mean profile to an appropriate value for the density bin.
bins = ed.sections.choose_bins(
natre.gamma_n.mean(["latitude", "longitude"]), np.arange(150, 2001, 10), decimals=6
)
props = natre[["temp", "salt", "gamma_n"]].interpolate_na("depth")
props["depth_"] = ("depth", props.depth.data)
mean_props = props.groupby_bins("gamma_n", bins).mean()
sliding = (
mean_props # .swap_dims({"gamma_n_bins": "depth_"})
# .rename({"depth_": "depth"})
.rolling(gamma_n_bins=10, center=True).construct("bin")
)
Δz = sliding.depth_.isel(bin=[-1, 0]).diff("bin").squeeze()
gradients = sliding.polyfit("bin", deg=1).sel(degree=1) / Δz
mean_props["Tz"] = gradients.temp_polyfit_coefficients
mean_props["Sz"] = gradients.salt_polyfit_coefficients
mean_props["N2"] = -9.81 / 1030 * gradients.gamma_n_polyfit_coefficients
mean_props
Mean gradients in a chosen O(100m) γ_n bin#
bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
binned = natre.reset_coords().groupby_bins("gamma_n", bins=bins)
i = 0
for label, group in binned:
i = i + 1
if i == 5:
break
ds = group.unstack()
ds
<xarray.Dataset>
Dimensions: (latitude: 10, longitude: 10, depth: 547)
Coordinates:
* latitude (latitude) float64 27.5 27.1 26.7 ... 24.7 24.3 23.9
* longitude (longitude) float64 -30.7 -30.3 -29.8 ... -27.2 -26.8
* depth (depth) float64 594.6 595.0 595.8 ... 495.8 497.8 499.8
Data variables: (12/18)
chi (latitude, longitude, depth) float64 1e-15 ... 6.867e-11
eps (latitude, longitude, depth) float64 7.443e-11 ... 2....
salt (latitude, longitude, depth) float64 35.56 ... 35.62
temp (latitude, longitude, depth) float64 11.7 11.7 ... 11.89
gamma_n (latitude, longitude, depth) float64 27.16 ... 27.16
pres (latitude, longitude, depth) float64 599.5 ... 503.5
... ...
chi_masked (latitude, longitude, depth) float64 1e-15 ... 6.867e-11
Krho (latitude, longitude, depth) float64 1.744e-06 ... 4....
KrhoTz (latitude, longitude, depth) float64 1.276e-08 ... 6....
eps_chi (latitude, longitude, depth) float64 3.988e-16 ... 1....
Kt (latitude, longitude, depth) float64 9.347e-12 ... 1....
KtTz (latitude, longitude, depth) float64 6.836e-14 ... 2....
Attributes: (12/13)
Conventions: CF-1.6
netcdf_version: 4
project: North Atlantic Tracer Release Experiment (NATRE)
expocode: 32OC250_4
cast_number: 3.0
title: Microstructure profiler data from the ship Oceanus...
... ...
latitude: 27.533166666666666
longitude: -30.723333333333333
chief_scientist: Raymond W. Schmitt
data_originator: Polzin
institution: WHOI
data_assembly_center: CCHDOed.sections.average_density_bin(group).dTdz_m
-0.011270900780997465
0.012996880869341362
<xarray.DataArray 'dTdz_m' ()>
array(0.0112709)
Attributes:
name: $∂_z θ_m$
units: °C/mgroup.unstack().depth
<xarray.DataArray 'depth' (depth: 547)> array([594.6, 595. , 595.8, ..., 495.8, 497.8, 499.8]) Coordinates: * depth (depth) float64 594.6 595.0 595.8 596.2 ... 493.4 495.8 497.8 499.8
ed.sections.fit1D(group, "theta", "depth", debug=True)
<xarray.DataArray 'depth' (depth: 10)>
array([551.590043, 561.407269, 572.258639, 583.53746 , 593.436723, 604.726001,
616.072442, 627.213435, 637.901493, 648.193224])
Coordinates:
gamma_n_bins (depth) object (27.16, 27.173] ... (27.277, 27.29]
* depth (depth) float64 551.6 561.4 572.3 583.5 ... 627.2 637.9 648.2
-0.011354461927299291
array(-0.01135446)
Computing along latitude and longitude lines#
def compute_single_line(section):
return ed.sections.bin_average_vertical(
section.cf.stack({"cast": ("latitude", "longitude")}).drop("cast"),
"neutral_density",
bins,
)
lonlines = natre.groupby("longitude", squeeze=False).map(compute_single_line)
latlines = natre.groupby("latitude", squeeze=False).map(compute_single_line)
latlines.to_netcdf("../datasets/natre-along-lat-lines.nc")
lonlines.to_netcdf("../datasets/natre-along-lon-lines.nc")
Buoyancy Reynolds Number#
There are a bunch of low values
natre["ν"] = dcpy.oceans.visc(natre.salt, natre.temp, natre.pres)
natre["κ"] = dcpy.oceans.tdif(natre.salt, natre.temp, natre.pres)
natre["Reb"] = natre.eps / natre.ν / natre.N2
natre["Reb"].attrs["long_name"] = "$Re_b$"
del natre.Reb.attrs["units"]
natre["χmol"] = 2 * natre.κ * natre.Tz**2
_, bins, _ = np.log10(natre.χmol + natre.chi).plot.hist(
bins=101, histtype="step", density=True, lw=2, aspect=3, size=3
)
np.log10(natre.chi).plot.hist(bins=bins, histtype="step", density=True, lw=2);
kwargs = dict(bins=np.linspace(-2, 5, 101), density=True, histtype="step", lw=2)
np.log10(natre.Reb.sel(latitude=24.7, method="nearest")).plot.hist(
**kwargs, aspect=3, size=3
)
np.log10(natre.Reb.sel(latitude=25.5, method="nearest")).plot.hist(**kwargs)
dcpy.plots.linex(np.log10([16, 200]))
/Users/dcherian/work/python/xarray/xarray/core/computation.py:734: RuntimeWarning: invalid value encountered in log10
result_data = func(*input_data)
/Users/dcherian/work/python/xarray/xarray/core/computation.py:734: RuntimeWarning: invalid value encountered in log10
result_data = func(*input_data)
NATRE \(ε\) vs \(ε_χ\)#
np.log10(
natre.eps_chi.where(np.abs(natre.Tz) > 5e-3).where(natre.eps_chi < 1e-6)
).plot.hist(bins=101, histtype="step", density=True)
np.log10(natre.eps).plot.hist(bins=101, histtype="step", density=True, yscale="linear");
natre.eps_chi.where(np.abs(natre.Tz) > 5e-3).where(natre.eps_chi < 5e-6).mean(
["latitude", "longitude"]
).coarsen(pres=100, boundary="trim").mean().plot()
natre.eps.mean(["latitude", "longitude"]).coarsen(
pres=100, boundary="trim"
).mean().plot(yscale="log")
[<matplotlib.lines.Line2D at 0x16f00e0a0>]
Compute NATRE microscale budget terms#
In theory, these bins are approx. 100m apart in neutral density
natre = ed.natre.read_natre(load=True).drop("depth")
bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
natre.gamma_n.stack({"latlon": ["latitude", "longitude"]}, create_index=False).cf.plot(
hue="latlon",
y="pres",
lw=0.5,
add_legend=False,
yincrease=False,
)
dcpy.plots.linex(bins)
micro.load(scheduler=client).to_netcdf("../datasets/natre-density-binned.nc")
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
T profiles#
ds = natre.stack(
{"station": ("latitude", "longitude")}, create_index=False
).assign_coords(station=np.arange(100))
import xfilter
Zname = ds.cf.axes["Z"][0]
CT = ds.CT.interpolate_na(Zname)
sortlo = xfilter.lowpass(CT, coord=Zname, freq=1 / 20)
sortlo100 = xfilter.lowpass(CT, coord=Zname, freq=1 / 100)
(
sortlo100.hvplot.line(y=Zname, color="k")
# * unsorted.CT.hvplot.line(y=Zname)
* sortlo.hvplot.line(y=Zname, color="r")
# * sort.CT.hvplot.line(y="pressure", ylim=(2000, 0))
).opts(ylim=(1200, 800), xlim=(6, 10), frame_width=500, frame_height=400)
Debugging estimates#
natre = ed.natre.read_natre(load=True, stack=False)
natre
13 40
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in get_gradient_sign
return self.ufunc(*args, **kwargs)
0 values == 0
28090
<xarray.Dataset>
Dimensions: (latitude: 10, longitude: 10, pres: 3981)
Coordinates:
* latitude (latitude) float64 27.5 27.1 26.7 26.3 ... 25.1 24.7 24.3 23.9
* longitude (longitude) float64 -30.7 -30.3 -29.8 -29.4 ... -27.6 -27.2 -26.8
* pres (pres) float64 10.0 10.5 11.0 11.5 ... 1.999e+03 2e+03 2e+03
time (latitude, longitude) datetime64[ns] 1992-03-28T15:28:59.99999...
depth (latitude, longitude, pres) float64 nan nan ... 1.978e+03
Data variables: (12/22)
chi (latitude, longitude, pres) float64 nan nan ... 1.145e-11
eps (latitude, longitude, pres) float64 nan nan ... 2.916e-11
salt (latitude, longitude, pres) float64 nan nan nan ... 35.05 35.05
temp (latitude, longitude, pres) float64 nan nan nan ... 3.934 3.933
gamma_n (latitude, longitude, pres) float64 nan nan nan ... 27.98 27.98
SA (latitude, longitude, pres) float64 nan nan nan ... 35.22 35.22
... ...
Kt (latitude, longitude, pres) float64 nan nan ... 2.579e-07
KtTz (latitude, longitude, pres) float64 nan nan ... 1.215e-09
Tz~ (latitude, longitude, pres) float64 nan nan nan ... nan nan nan
chi~ (latitude, longitude, pres) float64 nan nan nan ... nan nan nan
Kt~ (latitude, longitude, pres) float64 nan nan nan ... nan nan nan
KtTz~ (latitude, longitude, pres) float64 nan nan nan ... nan nan nan
Attributes: (12/13)
Conventions: CF-1.6
netcdf_version: 4
project: North Atlantic Tracer Release Experiment (NATRE)
expocode: 32OC250_4
cast_number: 3.0
title: Microstructure profiler data from the ship Oceanus...
... ...
latitude: 27.533166666666666
longitude: -30.723333333333333
chief_scientist: Raymond W. Schmitt
data_originator: Polzin
institution: WHOI
data_assembly_center: CCHDOfrom xarray.tests import raise_if_dask_computes
bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 150))
with raise_if_dask_computes():
micro150 = ed.sections.bin_average_vertical(
natre.stack(
{"cast": ["latitude", "longitude"]}, create_index=False
).assign_coords(cast=("cast", np.arange(100), {"cf_role": "profile_id"})),
"neutral_density",
bins,
blocksize=20,
)
2022-06-01 15:27:01,089 - distributed.nanny - WARNING - Restarting worker
2022-06-01 15:27:01,354 - distributed.nanny - WARNING - Restarting worker
micro150.load(client=client)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
from pandas import Int64Index as NumericIndex
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
from pandas import Int64Index as NumericIndex
<xarray.Dataset>
Dimensions: (bound: 2, gamma_n: 12, cast: 100)
Coordinates:
* bound (bound) <U5 'lower' 'upper'
* gamma_n (gamma_n) float64 26.61 26.91 27.12 ... 27.9 27.94 27.96
num_obs (gamma_n) int64 27405 28834 27525 ... 30729 29782 21047
pres (gamma_n) float64 225.0 375.1 ... 1.746e+03 1.878e+03
gamma_n_bounds (bound, gamma_n) float64 26.43 26.78 27.03 ... 27.95 27.97
latitude (cast) float64 27.5 27.5 27.5 27.5 ... 23.9 23.9 23.9 23.9
longitude (cast) float64 -30.7 -30.3 -29.8 ... -27.6 -27.2 -26.8
time (cast) datetime64[ns] 1992-03-28T15:28:59.999996928 ......
* cast (cast) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Data variables: (12/59)
chi (gamma_n) float64 1.538e-08 4.209e-09 ... 2.571e-10
eps (gamma_n) float64 9.01e-10 5.513e-10 ... 1.245e-10
KtTz (gamma_n) float64 1.978e-07 1.337e-07 ... 2.619e-08
KtTz~ (gamma_n) float64 2.604e-07 1.44e-07 ... 5.607e-08
chib2 (gamma_n) float64 7.689e-09 2.105e-09 ... 1.285e-10
hm (gamma_n) float64 146.5 152.9 146.5 ... 165.6 159.5 114.6
... ...
δKtTzTz (gamma_n) float64 1.479e-09 3.49e-10 ... 1.771e-11
δKtTz~Tz (gamma_n) float64 6.41e-10 3.147e-10 ... 1.952e-11
δwTTz (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
δresidual (gamma_n) float64 7.685e-10 2.61e-11 ... 9.237e-12
δresidual_chi (gamma_n) float64 2.518e-10 -4.472e-11 ... -2.887e-12
δpres (gamma_n) float64 73.26 76.45 73.26 ... 82.8 79.77 57.32
Attributes: (12/14)
Conventions: CF-1.6
netcdf_version: 4
project: North Atlantic Tracer Release Experiment (NATRE)
expocode: 32OC250_4
cast_number: 3.0
title: Microstructure profiler data from the ship Oceanus...
... ...
longitude: -30.723333333333333
chief_scientist: Raymond W. Schmitt
data_originator: Polzin
institution: WHOI
data_assembly_center: CCHDO
commit: eddydiff: b55f9526739dd63a768df2f816bc8327bcc0d247...Latest one with “sign stabilized” \(⟨K_T T_z⟩\) (purple). For some reason this one is sometimes larger than \(⟨χ⟩\) in the upper bins
ed.plot.debug_section_estimate(micro150, KρTz2=True)
micro150.to_netcdf("../datasets/natre-density-binned.nc")
Older#
Notation:
"KtTz"is “naive”"KtTz~"is estimated over 5m bins. I try to make the sign of the gradient somewhat “stable”."wTTz"is estimated in density space, without accounting for sign of \(\widetilde{T_z}\)
Latest one with “sign stabilized” \(⟨K_T T_z⟩\). For some reason this one is larger than \(⟨χ⟩\) in the upper bins
ed.plot.debug_section_estimate(micro150, KρTz2=True)
Bigger error bars with 100m bins. Averaging to 150m bins reduces the error.
dcpy.plots.fill_between_bounds(micro100, "KtTz", y="pres")
dcpy.plots.fill_between_bounds(micro150, "KtTz", y="pres")
plt.xscale("log")
dcpy.plots.fill_between_bounds(micro150, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro150, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro150, "KtTzTz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-12, 1e-8])
plt.xscale("log")
THe decrease in “naive” 0.5-m \(⟨K_T T_z⟩\) is because I was using \(T_z\) > 1e-3. This ends up throwing out a lot of points near the bottom.
dcpy.plots.fill_between_bounds(micro, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KtTzTz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-12, 1e-8])
plt.xscale("log")
The black is estimated in density space without accounting for sign of \(\widetilde{T_z}\), so we recover χ
dcpy.plots.fill_between_bounds(micro, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KtTzTz", y="pres")
dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-12, 1e-8])
plt.xscale("log")
Different bootstrapping method#
Explicitly calculate degrees of freedom like in Ferrari & Polzin (2005)
# micro.pres.attrs["bounds"] = "pres_err"
for label, group in grouped:
print(label)
break
(27.74, 27.79]
ds = group.unstack()
da = ds.KtTz
dp = 0.5
corrscales = [[slice(0, 800), 5], [slice(800, None), 10]]
dof = ed.sections.get_dof(da, dp=0.5, corrscales=corrscales)
dof
1011
from arch.bootstrap import IIDBootstrap
array = da.data.reshape(-1)
print("old...", ed.sections.compute_bootstrapped_mean_ci(array, 20, clean=False))
array = array[~np.isnan(array)]
plt.plot(array, marker="x")
absarray = array # np.abs(array)
thresh = np.mean(absarray) + 30 * np.std(absarray)
array = np.where(np.abs(array) < thresh, array, np.nan)
array = array[~np.isnan(array)]
plt.plot(array, marker="x")
mean = np.mean(array)
stderr = np.mean(IIDBootstrap(array[~np.isnan(array)]).apply(np.std)) / np.sqrt(dof)
print("new...", [mean - stderr, mean, mean + stderr])
ed.sections.compute_bootstrapped_mean_ci(array, 20, clean=False)
old... [4.50550872e-08 6.73010347e-08 1.42173197e-07]
new... [1.784538731438142e-08, 4.454195776804866e-08, 7.12385282217159e-08]
array([3.21298005e-08, 4.45419578e-08, 6.02803924e-08])
Partitioning \(θ_z\)#
def bin_coarsen(ds, func="mean"):
return getattr(
ds.coarsen(latitude=10, longitude=10, pres=300, boundary="trim"), func
)()
natre = ed.natre.read_natre(load=True).sel(pres=slice(2000))
natre_binned = xr.load_dataset("../datasets/natre-density-binned.nc")
natre_coarse = natre.where(natre.chi_masked.notnull()).pipe(bin_coarsen).squeeze()
f, ax = plt.subplots(1, 4, sharey=False)
for Tzmin in [1e-3, 2e-3, 5e-3]:
mask = (np.abs(natre.Tz) > Tzmin) & (natre.chi < 1e-6)
chib2 = (natre.chi / 2).where(mask)
Kt = chib2 / natre.Tz.where(mask) ** 2
KtTz = chib2 / natre.Tz.where(mask)
bin_coarsen(chib2).squeeze().cf.plot(
y="pres", label=str(Tzmin), ax=ax[0], xscale="log"
)
bin_coarsen(Kt).squeeze().cf.plot(
y="pres", label=str(Tzmin), ax=ax[1], xscale="log"
)
bin_coarsen(KtTz).squeeze().cf.plot(
y="pres", label=str(Tzmin), ax=ax[2], xscale="log"
)
(bin_coarsen(KtTz, func="std") / np.sqrt(100 * 10)).squeeze().cf.plot(
y="pres", label=str(Tzmin), ax=ax[3], xscale="log"
)
# np.log(1025 * 4000 * np.abs(KtTz)).plot.hist(
# bins=np.arange(-10, 10, 0.05), histtype="step", ax=ax[3], density=True
# )
(natre_binned.chib2).cf.plot(y="pres", ax=ax[0], color="k")
(natre_binned.KtTzTz / natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[2])
(natre_binned.chib2 / natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[2], color="k")
(natre_binned.Krho_m * natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[2], color="c")
(natre_binned.δKtTzTz / natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[3], color="k")
(natre_binned.δKρTz2 / natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[3], color="k")
ax[1].legend()
dcpy.plots.clean_axes(ax)
f.set_size_inches((9, 4))
Kt = natre.chi / 2 / natre.Tz.where(np.abs(natre.Tz) > 5e-3) ** 2
natre_coarse.Kt.plot(y="pres")
(natre_binned.Krho_m).plot(y="pres", color="k")
(natre_binned.Kt_m / natre_binned.dTdz_m).cf.plot(y="pres", xscale="log")
bin_coarsen(Kt).plot(y="pres")
# Kt.coarsen({"latitude": 10, "longitude": 10, "pres_3.1": 30}, boundary="trim").mean().plot(y="pres_3.1")
[<matplotlib.lines.Line2D at 0x16920d910>]
(bin_coarsen(Kt) * natre_coarse.temp.differentiate("pres") ** 2).squeeze().cf.plot()
(natre_coarse.Kt * natre_coarse.temp.differentiate("pres") ** 2).cf.plot(xscale="log")
(natre_coarse.chi / 2).cf.plot(xscale="log")
natre_binned.chib2.plot(y="pres")
[<matplotlib.lines.Line2D at 0x168e38100>]
dTdz.coarsen(
{"pres_3.1": 61, "latitude": 10, "longitude": 10}, boundary="trim"
).mean().plot()
natre_binned.dTdz_m.plot(x="pres", yscale="log")
[<matplotlib.lines.Line2D at 0x178c62370>]
natre_binned.dTdz_m.plot(y="pres")
(-1 * natre_coarse.temp.differentiate("pres")).plot(y="pres")
[<matplotlib.lines.Line2D at 0x16a09be20>]
Compare neutral density bins to FP2005#
plt.figure()
fp = xr.DataArray(
[26.59, 26.95, 27.22, 27.46, 27.64, 27.77, 27.85, 27.91, 27.95],
dims=("depth",),
coords={"depth": [200, 400, 600, 800, 1000, 1200, 1400, 1600, 1800]},
)
fp.plot(marker="x", label="FP2005")
micro.gamma_n.plot(x="pres", marker=".", label="mine")
# plt.plot(np.arange(150, 2001, 100), bins, marker='^')
plt.gca().set_xticks(fp.depth)
plt.legend()
<matplotlib.legend.Legend at 0x7f2cd6b69f70>
Better \(\widetilde{K_T T_z}\) estimate#
It appears that preserving the sign of \(wT\) is important. After a lot of averaging I do get nearly all positive values but the error bars are big!
Yeah the sign for \(T_m + T_e\) can be negative and opposite to \(T_z^m\) Indeed that is important!
If I sort the \(Θ\) then I get rid of the stable inversions that are important.
~No matter what I do, I basically recover χ :(.~ I recover χ if I don’t pay attention to this fact.
I think this explains why \(\widetilde{K_t T_z}\) with sign preserved is closer to \(Γ ⟨ε⟩/N_m^2 T_z^m\). The cancellations are important! Using absolute values kills the cancellations…
Might be better to stabilize the sign somehow; Logic here is that
the sign of \(\widetilde{T_z}\) is set by the intermediate scale (30-50m-ish);
but we want \(\widetilde{w_tθ_t} = \widetilde{χ}/\widetilde{T_z}/2\) over 5-m scales
So filter at 10m wavelength; choose sign;
Then determine \(\widetilde{T_z}\) magnitude over 5-m linear fits;
apply sign of first element in window (this is somewhat arbitrary)
natre = ed.natre.read_natre(load=True)
natre_binned = xr.load_dataset("../datasets/natre-density-binned.nc")
Tzmi = natre_binned.dTdz_m.reset_coords(drop=True).interp(gamma_n=natre.gamma_n)
Tmi = natre_binned.theta.reset_coords(drop=True).interp(gamma_n=natre.gamma_n)
Depth-space#
Here I estimate \(\widetilde{w^tθ^t} = \widetilde{χ} / 2 / \widetilde{T}_z\) over 5m scales.
Using linear fit to \(T\) (or sorted \(T\)) ofr gradient;
simple 5m average of χ
The core problem is that a simple averaging of \(\widetilde{T}_z\) basically yields \(T_z^m\); but it should be bigger for this strategy to work.
(natre.Tz - Tzmi).plot.line(row="latitude", col="longitude")
<xarray.plot.facetgrid.FacetGrid at 0x17cf48d60>
import hvplot.xarray
import xfilter
loc = dict(latitude=3, longitude=8)
T = natre.CT.isel(loc).interpolate_na("pres").dropna("pres")
# Tm = xfilter.lowpass(T, coord="pres", freq=1 / 200, order=1)
Te = xfilter.lowpass(T, coord="pres", freq=1 / 20, order=1)
# (T.plot())
# (Te.plot())
# (Tmi.isel(loc).plot())
# xfilter.lowpass(
# (T - Tmi.isel(loc)).interpolate_na("pres"), coord="pres", freq=1 / 20
# ).plot()
# (Tm.plot())
# ((T - Te).plot(ax=plt.gca().twinx()))
# (T.hvplot.line(x="pres") * Te.hvplot.line(x="pres"))
(T.differentiate("pres").hvplot.line()) * Te.differentiate("pres").hvplot.line()
clean.dTdz.isel(latitude=2, longitude=6).plot()
clean.dTdz.isel(latitude=9, longitude=8).plot()
natre_binned.dTdz_m.plot(x="pres", color="k")
[<matplotlib.lines.Line2D at 0x17c902df0>]
(clean.dTdz).coarsen(
{"latitude": 10, "longitude": 10, "pres_5.1": 20}, boundary="trim"
).mean().plot(yscale="log")
(natre_binned.dTdz_m).plot(x="pres")
[<matplotlib.lines.Line2D at 0x17757ff40>]
clean = ed.sections.estimate_microscale_stirring_depth_space(
natre, filter_len=10, segment_len=3
)
plt.figure()
clean["Tz~"].coarsen(pres=200, latitude=10, longitude=10, boundary="trim").mean().plot()
natre_binned.dTdz_m.plot(x="pres", color="k", yscale="log")
plt.figure()
natre.KtTz.where(np.abs(natre.Tz) > 1e-5).mean(["latitude", "longitude"]).plot(
yscale="log", lw=0.5, label="simple"
)
(
(
natre.KtTz.where(np.abs(natre.Tz) > 1e-5)
.coarsen(pres=200, latitude=10, longitude=10, boundary="trim")
.mean()
).plot(yscale="log", label="simple coarsened")
)
(
clean["KtTz~"]
# .where(np.abs(clean["KtTz~"]) < 5e-6)
.cf.coarsen({"Z": 200, "latitude": 10, "longitude": 10}, boundary="trim").mean()
# .mean(["latitude", "longitude"])
.plot(label="cleaner")
)
(natre_binned.KρTz2 / natre_binned.dTdz_m).plot(label="$K_ρT_z^m$", color="k", x="pres")
(natre_binned.chib2 / natre_binned.dTdz_m).plot(
label="$χ/2/T_z^m$",
color="k",
x="pres",
ylim=(1e-8, 1e-6),
)
# (micro.wTTz / micro.dTdz_m).plot(label="$K_tT_z$", color="cyan", x="pres")
plt.legend()
7 20
50243 values < 0
<matplotlib.legend.Legend at 0x177280df0>
kwargs = dict(pres=100, latitude=10, longitude=10, boundary="trim")
clean.KtTz.mean(["latitude", "longitude"]).coarsen(
{"pres_5.1": 20}, boundary="trim"
).mean().plot()
(
((natre.KtTz).where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean()).plot(
yscale="log"
)
)
((natre.KrhoTz).coarsen(**kwargs).mean().plot(yscale="log", color="k"))
# np.abs(natre.chi / 2).mean(["latitude", "longitude"]).coarsen(**kwargs).mean().plot(
# yscale="log"
# )
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [316], in <cell line: 3>()
1 kwargs = dict(pres=100, latitude=10, longitude=10, boundary="trim")
----> 3 clean.KtTz.mean(["latitude", "longitude"]).coarsen(
4 {"pres_5.1": 20}, boundary="trim"
5 ).mean().plot()
7 (
8 (((natre.KtTz).where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean())).plot(
9 yscale="log"
10 )
11 )
13 ((natre.KrhoTz).coarsen(**kwargs).mean().plot(yscale="log", color="k"))
NameError: name 'clean' is not defined
natre.eps.coarsen(**kwargs).mean().plot(yscale="log")
natre.eps_chi.where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean().plot(yscale="log")
[<matplotlib.lines.Line2D at 0x16469f1f0>]
Interestingly excluding small gradients excludes small χ, so the \(\widetilde{K_TT_z}\) estimate increases with more filtering!
kwargs = dict(bins=100, density=True, histtype="step")
np.log10(chi / 2 / dTdz).plot.hist(**kwargs)
np.log10(chi / 2 / dTdz.where(dTdz > 5e-3)).plot.hist(**kwargs);
T-space#
This is not as easy.
Choosing T bins is a real pain. Because we are covering a giant region; the spread in cast-mean T is large. So picking bins is a little challening.
By averaging over all casts I end up with \(⟨χ⟩\), and without accounting for sign of \(Δz\) (separation between isotherms), I end up with \(T_z^m\) in the denomoinator.
natre = ed.natre.read_natre(load=True)
clean = ed.sections.estimate_microscale_stirring_depth_space(
natre.stack({"cast": ("latitude", "longitude")}).drop("cast"), # .isel(cast=36),
filter_len=10,
segment_len=5,
)
clean
11 20
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in get_gradient_sign
return self.ufunc(*args, **kwargs)
0 values == 0
33882
<xarray.Dataset>
Dimensions: (pres: 3981, cast: 100)
Coordinates:
* pres (pres) float64 10.0 10.5 11.0 11.5 ... 1.999e+03 2e+03 2e+03
* cast (cast) int64 0 1 2 3 4 5 6 7 8 9 ... 90 91 92 93 94 95 96 97 98 99
window (pres) float64 -0.0 -0.5 -1.0 -1.5 -2.0 ... nan nan nan nan nan
Data variables:
Tz~ (pres, cast) float64 nan nan nan nan nan ... nan nan nan nan nan
chi~ (pres, cast) float64 nan nan nan nan nan ... nan nan nan nan nan
Kt~ (pres, cast) float64 nan nan nan nan nan ... nan nan nan nan nan
KtTz~ (pres, cast) float64 nan nan nan nan nan ... nan nan nan nan nan
T~ (cast, pres) float64 nan nan nan nan nan ... nan nan nan nan nannatre = natre.update(clean)
from xarray.tests import raise_if_dask_computes
bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
with raise_if_dask_computes():
micro = ed.sections.bin_average_vertical(
natre.cf.stack(
{"cast": ("latitude", "longitude")}, create_index=False
).assign_coords({"cast": ("cast", np.arange(100), {"cf_role": "profile_id"})}),
"neutral_density",
bins,
blocksize=10,
)
# micro.pres.attrs["bounds"] = "pres_err"
micro
<xarray.Dataset>
Dimensions: (bound: 2, gamma_n: 18, cast: 100)
Coordinates:
* bound (bound) <U5 'lower' 'upper'
* gamma_n (gamma_n) float64 26.56 26.78 26.95 ... 27.93 27.95 27.96
num_obs (gamma_n) int64 19358 19463 20127 ... 20191 11530 22215
pres (gamma_n) float64 201.4 300.6 401.1 ... 1.795e+03 1.877e+03
gamma_n_bounds (bound, gamma_n) float64 26.43 26.69 26.87 ... 27.95 27.97
latitude (cast) float64 27.5 27.5 27.5 27.5 ... 23.9 23.9 23.9 23.9
longitude (cast) float64 -30.7 -30.3 -29.8 -29.4 ... -27.6 -27.2 -26.8
time (cast) datetime64[ns] 1992-03-28T15:28:59.999996928 ... 1...
* cast (cast) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Data variables: (12/54)
chi (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
eps (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
KtTz (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
KtTz~ (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
hm (gamma_n) float64 98.77 98.73 101.8 ... 101.3 58.27 114.8
wT (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
... ...
δKtTz~Tz (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
δwTTz (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
δresidual (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
δpres (gamma_n) float64 49.38 49.37 50.92 ... 50.63 29.14 57.4
chib2 (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
chib2_err (gamma_n, bound) float64 dask.array<chunksize=(1, 2), meta=np.ndarray>
Attributes: (12/14)
Conventions: CF-1.6
netcdf_version: 4
project: North Atlantic Tracer Release Experiment (NATRE)
expocode: 32OC250_4
cast_number: 3.0
title: Microstructure profiler data from the ship Oceanus...
... ...
longitude: -30.723333333333333
chief_scientist: Raymond W. Schmitt
data_originator: Polzin
institution: WHOI
data_assembly_center: CCHDO
commit: eddydiff: 3c5e37fac5bc4c53f9567addc4c058ef415a63b4...micro.load()
<xarray.Dataset>
Dimensions: (bound: 2, gamma_n: 18, cast: 100)
Coordinates:
* bound (bound) <U5 'lower' 'upper'
* gamma_n (gamma_n) float64 26.56 26.78 26.95 ... 27.93 27.95 27.96
num_obs (gamma_n) int64 19358 19463 20127 ... 20191 11530 22215
pres (gamma_n) float64 201.4 300.6 401.1 ... 1.795e+03 1.877e+03
gamma_n_bounds (bound, gamma_n) float64 26.43 26.69 26.87 ... 27.95 27.97
latitude (cast) float64 27.5 27.5 27.5 27.5 ... 23.9 23.9 23.9 23.9
longitude (cast) float64 -30.7 -30.3 -29.8 -29.4 ... -27.6 -27.2 -26.8
time (cast) datetime64[ns] 1992-03-28T15:28:59.999996928 ... 1...
* cast (cast) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Data variables: (12/54)
chi (gamma_n) float64 1.664e-08 5.465e-09 ... 2.303e-10
eps (gamma_n) float64 9.838e-10 5.747e-10 ... 1.235e-10
KtTz (gamma_n) float64 2.515e-07 1.472e-07 ... 3.07e-08 2.191e-08
KtTz~ (gamma_n) float64 2.887e-07 1.79e-07 ... 5.764e-08 5.291e-08
hm (gamma_n) float64 98.77 98.73 101.8 ... 101.3 58.27 114.8
wT (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
... ...
δKtTz~Tz (gamma_n) float64 7.741e-10 4.478e-10 ... 2.664e-11 1.7e-11
δwTTz (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
δresidual (gamma_n) float64 7.387e-10 5.777e-11 ... 6.443e-12
δpres (gamma_n) float64 49.38 49.37 50.92 ... 50.63 29.14 57.4
chib2 (gamma_n) float64 8.32e-09 2.733e-09 ... 1.428e-10 1.151e-10
chib2_err (gamma_n, bound) float64 7.713e-09 8.926e-09 ... 1.229e-10
Attributes: (12/14)
Conventions: CF-1.6
netcdf_version: 4
project: North Atlantic Tracer Release Experiment (NATRE)
expocode: 32OC250_4
cast_number: 3.0
title: Microstructure profiler data from the ship Oceanus...
... ...
longitude: -30.723333333333333
chief_scientist: Raymond W. Schmitt
data_originator: Polzin
institution: WHOI
data_assembly_center: CCHDO
commit: eddydiff: 3c5e37fac5bc4c53f9567addc4c058ef415a63b4...Using fancy sign detection#
def notnull(data):
flat = data.data.ravel()
return flat[~np.isnan(flat)]
mask = np.abs(clean["KtTz~"]) > 1e-5
clean["chi~"].where(mask).pipe(notnull)
array([2.29729636e-05, 4.99936064e-06, 8.38180596e-08, 4.42592055e-08,
2.45233000e-08, 6.14768901e-08, 5.21607218e-07, 2.45605935e-07,
9.83617336e-08, 2.11507639e-07, 1.67454668e-08])
clean["Tz~"].where(np.abs(clean["KtTz~"]) > 1e-5).pipe(notnull)
array([ 0.14074649, 0.13337417, 0.00319876, 0.00175999, 0.00105803,
0.00300432, 0.00974087, 0.00935505, -0.00478368, 0.00899659,
0.00078514])
clean["Tz~"].isel(cast=16).plot(marker="x")
[<matplotlib.lines.Line2D at 0x191bccb80>]
clean["KtTz~"].isel(cast=26).ffill("pres").plot(yscale="linear")
[<matplotlib.lines.Line2D at 0x193737850>]
ed.plot.debug_section_estimate(micro)
dcpy.plots.fill_between_bounds(micro, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KtTz~Tz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "KtTzTz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-11, 1e-8])
plt.xscale("log")
Single bin#
from xarray.tests import raise_if_dask_computes
bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
with raise_if_dask_computes():
grouped = ed.sections.bin_average_vertical(
natre.cf.stack(
{"cast": ("latitude", "longitude")}, create_index=False
).assign_coords({"cast": ("cast", np.arange(100), {"cf_role": "profile_id"})}),
"neutral_density",
bins[10:],
blocksize=10,
return_group=True,
)
# micro.pres.attrs["bounds"] = "pres_err"
for label, group in grouped:
break
Tbins, sortT = ed.sections.estimate_microscale_stirring(
group.unstack(), dz=7, debug=True
)
dT = np.diff(Tbins).mean()
sort = dcpy.oceans.thorpesort(sortT, by="gamma_n", core_dim="pres")
sort.CT.plot.line(hue="cast", add_legend=False);
sort.CT.isel(cast=[20, 30, 60]).cf.dropna("Z").plot(hue="cast")
dcpy.plots.liney(Tbins)
CT = profile.CT.data
profiles = natre.copy(deep=True)
profiles["CT"] = profiles.CT.interpolate_na("pres")
profiles["Tzfilt"] = -1 * xfilter.lowpass(
profiles.CT, coord="pres", freq=1 / 10
).differentiate("pres").bfill("pres")
stacked = profiles.stack({"cast": ("latitude", "longitude")}, create_index=False)
profile = stacked.interpolate_na("pres") # .sel(pres=slice(1120, 1250))
# get_gradient_sign(
# profile.pres.data, profile.CT.data, profile.Tzfilt.data, 5, 10, a.data
# )
# profile.CT.plot()
Tzsign.plot(ax=plt.gca().twinx())
<matplotlib.collections.QuadMesh at 0x185d14790>
def clean_idx(idx, CT, ΔP, window):
P = CT["pres"].data
CT = CT.data
def _clean(idx, P, ΔP, window):
newidx = [idx[0]]
for left, right in zip(idx[:-1], idx[1:]):
Pl, Pr = P[left], P[right]
# print(np.abs(Pr - Pl))
if np.abs(Pr - Pl) > ΔP:
if ((right - left) % window) < window // 2:
right += 1
newidx.append(right)
if newidx[-1] != idx[-1]:
# Last index was dropped
# instead to preserve all data;
# we drop the previous one and re-add the last one
newidx.pop()
newidx.append(idx[-1])
print(idx, newidx)
return newidx
idx = _clean(idx, P, ΔP, window=window)
idx = _clean(idx, CT, 1e-4 * ΔP, window=1e10)
return idx
def smooth_gradient(profile, chi, ΔP, window, debug=False):
dp = 0.5
smooth = xfilter.lowpass(
profile, coord="pres", freq=1 / int(ΔP // dp), num_discard=5
)
Tzsmooth = smooth.differentiate("pres")
# array = profile.data
if debug:
f, ax = plt.subplots(4, 1, sharex=True, constrained_layout=True)
profile.plot(ax=ax[0])
smooth.plot(ax=ax[0])
(idx,) = np.nonzero(np.abs(np.sign(Tzsmooth).diff("pres").data) > 0)
idx = np.concatenate([[0], idx, [profile.sizes["pres"] - 1]])
if debug:
dcpy.plots.linex(Tzsmooth["pres"][idx], color="k", lw=0.5, ax=ax[0])
idx = clean_idx(idx, profile, ΔP=ΔP, window=window)
if debug:
dcpy.plots.linex(Tzsmooth["pres"][idx], color="b", ax=ax[0])
Tzs = []
for subprof in np.split(profile, idx[1:-1]):
# -1 for temperature slope
sign = np.sign(subprof[0] - subprof[-1]).data
coarse = subprof.coarsen(pres=window, boundary="pad").construct(
pres=("p", "window")
)
sorty = dcpy.oceans.thorpesort(
coarse, coarse, core_dim="window", ascending=False
)
sorty = sorty.where(sorty.count("window") > 2)
Tzsubprof = sign * np.abs(
sorty.polyfit("window", deg=1)
.sel(degree=1, drop=True)
.polyfit_coefficients.data
)
Tzs.append(np.repeat(Tzsubprof, window)[..., : subprof.shape[-1]])
dTdz = profile.copy(data=np.concatenate(Tzs)).ffill("pres")
if debug:
dTdz.plot(ax=ax[1])
(-1 * Tzsmooth).plot(ax=ax[1])
(chi).plot(yscale="log", ylim=(1e-12, None), ax=ax[2])
ax2 = ax[2].twinx()
(dTdz).plot(ax=ax2, color="k")
(dTdz.where(np.abs(dTdz) < 2e-4)).plot(color="b", marker="o", ls="none", ax=ax2)
dcpy.plots.set_axes_color(ax2, "k")
(1025 * 4000 * -chi / dTdz).plot(ax=ax[3], yscale="symlog")
(1025 * 4000 * -chi / -Tzsmooth).plot(ax=ax[3])
(1025 * 4000 * -chi / dTdz.where(np.abs(dTdz) < 2e-4)).plot(
color="b", marker="o", ls="none", ax=ax[3]
)
dcpy.plots.clean_axes(np.atleast_2d(ax).T)
f.set_size_inches((10, 7))
return dTdz
# dcpy.plots.liney(Tbins)
# NATRE 36, 56, 60 are good to test with
profile = sort.isel(cast=36).cf.dropna("Z")
dTdz = smooth_gradient(profile.CT, profile.chi, ΔP=5, window=20, debug=True);
[ 0 6 11 37 59 76 104 119 124 135 148 153 187] [0, 38, 60, 76, 105, 119, 135, 148, 187]
[0, 38, 60, 76, 105, 119, 135, 148, 187] [0, 39, 61, 77, 106, 120, 136, 149, 187]
Tz_m = sort.CT.polyfit("pres", deg=1).sel(degree=1, drop=True).polyfit_coefficients
Tz_m.plot()
[<matplotlib.lines.Line2D at 0x180d0ceb0>]
Tz = -1 * xfilter.lowpass(sort.CT, coord="pres", freq=1 / 10).differentiate("pres")
KtTz = sort.chi / Tz
KtTz = KtTz.where(np.abs(KtTz) < 300 / 1025 / 4000)
(KtTz.mean("pres") * Tz_m).plot()
(sort.chi.mean("pres") / 2).plot()
[<matplotlib.lines.Line2D at 0x1814e42b0>]
Try ε#
natre.eps_chi.where(np.abs(natre.Tz) > 1e-10).mean(["latitude", "longitude"]).plot(
yscale="log"
)
natre.eps_chi.where(np.abs(natre.Tz) > 1e-2).mean(["latitude", "longitude"]).plot(
yscale="log"
)
natre.eps.mean(["latitude", "longitude"]).plot(yscale="log")
[<matplotlib.lines.Line2D at 0x170f51220>]
Sensitivity of χpod estimate to mean K estimate#
Various estimates of the variance produced by turbulent stirring of mean vertical gradient
Ferrai & Polzin (2005): \(Γ ⟨ε⟩/N_m^2 ∂_zθ_m^2\)
χpod estimate: \(⟨K_T θ_z⟩ ∂_zθ_m\)
(1) with \(ε_χ\) : \(Γ ⟨ε_χ⟩/N_m^2 ∂_zθ_m^2\)
(2) with \(K_ρ\) instead of \(K_T\): \(⟨K_ρ θ_z⟩ ∂_zθ_m^2\)
(chidens.Krho_m * chidens.dTdz_m**2).cf.plot.step(
color="k", lw=2, label="$Γ ⟨ε⟩/N_m^2 ∂_zθ_m^2$"
)
(chidens.KtTz * chidens.dTdz_m).cf.plot.step(
color="k", marker="x", ls="none", label="$⟨K_T θ_z⟩ ∂_zθ_m$"
)
(0.2 * chidens.eps_chi / chidens.N2_m * chidens.dTdz_m**2).cf.plot.step(
color="C1", lw=2, label="$Γ ⟨ε_χ⟩/N_m^2 ∂_zθ_m^2$"
)
(chidens.KrhoTz * chidens.dTdz_m).cf.plot.step(
lw=2, label="$⟨K_ρ θ_z⟩ ∂_zθ_m$", color="C2"
)
plt.grid(True, which="both", lw=0.5)
plt.legend()
plt.xscale("log")
plt.xlabel("Variance production or dissipation [°C²/s]")
plt.gcf().set_size_inches((4, 5))
f, ax = plt.subplots(1, 2, sharey=True, sharex=True)
micro.Krho.cf.plot.step(y="pres", xscale="log", ax=ax[0])
micro.Kt.cf.plot.step(y="pres", ax=ax[0])
micro.Krho_m.cf.plot.step(y="pres", ax=ax[1])
micro.Kt_m.cf.plot.step(y="pres", ax=ax[1])
ax[0].legend(["$⟨K_ρ⟩$", "$⟨K_T⟩$"])
ax[1].legend(["$K_ρ^m = Γ ⟨ε⟩/N_m²$", "$K_T^m = ⟨χ⟩/2θ_m^2$"])
[aa.grid(True, which="both", lw=0.5) for aa in ax]
plt.gcf().set_size_inches((6, 5))